Physiological and transcriptome analyses highlight multiple pathways involved in drought stress in Medicago falcata

Medicago falcata is one of the leguminous forage crops, which grows well in arid and semiarid region. To fully investigate the mechanism of drought resistance response in M. falcata, we challenged the M. falcata plants with 30% PEG-6000, and performed physiological and transcriptome analyses. It was found that, the activities of antioxidant enzymes (eg. SOD, POD, and CAT) and soluble sugar content were all increased in the PEG-treated group, as compared to the control group. Transcriptome results showed that a total of 706 genes were differentially expressed in the PEG-treated plants in comparison with the control. Gene enrichment analyses on differentially expressed genes revealed that a number of genes in various pathway were significantly enriched, including the phenylpropanoid biosynthesis (ko00940) and glycolysis/gluconeogenesis (ko00010), indicating the involvement of these key pathways in drought response. Furthermore, the expression levels of seven differentially expressed genes were verified to be involved in drought response in M. falcata by qPCR. Taken together, these results will provide valuable information related to drought response in M. falcata and lay a foundation for molecular studies and genetic breeding of legume crops in future research.


Introduction
Water shortage is one of the most serious abiotic factors affecting plant growth and production, and global climate change threatens sustainable agriculture and may increase the frequency and severity of drought, which in turn prompt to develop drought-tolerant varieties of many crops [1][2][3][4]. Plants are very vulnerable to drought stress, and water shortage is fatal to plants and leads to huge economic losses. During the long-term evolution, plants have acquired a series of response protection mechanisms to endure drought stress [5], including morphological, physiological, biochemical and molecular responses. At cellular level, drought signals promote production of stress-protectant metabolites (eg. proline and soluble sugar), drought resistance of M. falcata by RNA-seq is the most critical endeavors to develop necessary tools and integrate molecular breeding approaches in alfalfa. In this study, we explored transcriptome sequencing technique to enrich the transcriptome in M. falcata and identified genes related to drought resistance. We also screened a number of differentially expressed genes related drought resistance and plant recovery. The physiological responses and transcriptome data of drought resistance can provide information for the breeding of drought resistance in both M. falcata and alfalfa.

Plant materials and drought treatment
M. falcata seeds were collected in Xinjiang, and stored at the Grassland Resource and Ecology Key Laboratory of Xinjiang Agricultural University. After physical scarification, seeds were placed on petri dishes and incubated at 4˚C for 24h before germinated at 25˚C (16h light/8h dark). After germination, seedlings were grown in 3:1:1 topsoil/vermiculite/perlite (v/v/v) mixture at 25˚C/18˚C (16h/8h day/night). After growth for 1 month, 30% PEG-6000 was added to simulate drought stress, and water was used as control. After one hour PEG treatment, the seedlings (about 60 plants) were harvested as treatment group, and seedlings with water were harvested as control group. The remaining seedlings (about 30) were then subjected to a rehydration treatment for one hour as recovery group. Each group contained a triplicate with more than 3 seedlings for each replicate, 9 samples in total. The whole seedlings (about 2.0g) were immediately rinsed in water, dried with filter paper, and timely placed in liquid nitrogen, and stored at -80˚C for RNA extraction and transcriptome analysis (CK1, CK2, CK3 for the control group; T1, T2, T3 for the PEG-treatment group; R1, R2, R3 for the recovery group). Simultaneously, the roots and leaves were further separated, frozen in liquid nitrogen, and stored at -80˚C for physiological index measurement, three biological replicates were set for each control group, PEG treatment group and recovery group.

Measurement of physiological indicators
Half grams of the leaves and root samples were used to quantify relative conductivity as described by Ding and Wang [29], and soluble sugar content by anthrone colorimetry [30]. The activity of antioxidant enzymes including superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT) were assayed using the nitrogen blue tetrazolium (NBT) method [31], the guaiacol method [32], and the UV absorption method [33], respectively. Data were analyzed using one-way ANOVA with software SPSS 20.0.

Sequence analysis, assembly, and annotation
Raw data (reads) in fastq format were firstly processed using an in-house Perl scripts. In this step, clean reads were obtained by removing reads containing adapter sequences, poly-N sequences, or low-quality reads. Then the Q20, Q30, and GC content and sequence duplication were determined, and de novo assembly of the transcriptome was carried out using Trinity, with min_kmer_cov value of 2 and all the other parameters were default [34]. Assembly quality was assessed by the length distribution of unigenes. The read counts of unigenes were calculated by the software RSEM v1.2.15 [35], and the gene expression level was calculated using expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced (FPKM). Plant transcription factor genes were predicted by using iTAK 1.2 software. The basic principle is to use the well-defined transcription factor family genes in the database to identify TF in M. falcata through hmmscan. Using RSEM v1.2.15 software, unigenes were mapped with the reference generated by Trinity software, functions of unigenes were annotated by BLAST queries of the NCBI non-redundant protein (Nr) database, the NCBI non-redundant nucleotide sequences (Nt) database, the protein families (Pfam) database, the Swiss-Prot protein database, the Kyoto Encyclopedia of Genes and Genomes (KEGG) Ortholog database, the Gene Ontology (GO) database, and the Clusters of Orthologous Groups of proteins (KOG) database. All queries were performed using a cut-off E-value of 1×10 −5 .

Gene expression and enrichment analysis
Analysis on differential expression genes of the treatment and control groups was performed using the DEseq R package (1.10.1), and applied with read count values to analyze gene expression levels with a threshold of |log 2 Fold change|�1 and padj�0.05. The DEGs were visualized with Venn diagrams and Volcano plots. Clustering was performed using the pheatmap R package. GO enrichment analysis of DEGs was carried out using Goseq [36] with corrected P<0.05 as threshold. KEGG enrichment analysis of DEGs was performed using KOBAS version 2.0, BH correction for FDR. KEGG pathways with q value�0.05 were considered to be significantly enriched DEGs, and the top 20 pathway entries were presented with enrichment scatter plot.

qPCR analyses
Quantitative real-time PCR (qPCR) analysis was performed to validate the RNA-Seq results. The primers (S1 Table) were designed for the selected 10 genes using Primer Premier 5.0 Tool, and actin gene was used as internal control [37]. Reverse transcription was performed with total RNA using M-MLV Reverse Transcriptase (Promega, Madison, WI, USA). qPCR analysis was conducted using a CFX-96 Real-Time System (Bio-Rad, Hercules, CA, USA), and SYBR Premix Ex Taq (TaKaRa, Kyoto, Japan). qPCRs were performed using the following conditions: denaturation at 95˚C for 15 min, followed by 40 cycles of amplification (95˚C for 10s, 58˚C for 20s, and 72˚C for 20 s). Each sample was performed with a technical triplicates, and gene expression was evaluated using the 2− ΔΔCt method.

Physiological changes of M. falcata under drought treatment
We used PEG solution to mimic drought stress, and one-month-old M. falcata seedlings were treated with PEG solution, with seedlings treated with water as control. Various physiology index from both leaves and roots samples were analyzed from treatment and control seedlings. It was revealed that, under drought stress, conductivity and soluble sugar content were increased significantly (P<0.05), whether in leaves or in roots (Fig 1A and 1B). Meanwhile, the enzyme activity of SOD, POD, and CAT were all increased at various level in comparison with the untreated control in both leaves and roots (Fig 1C-1E). After re-watering, seedlings of drought treatment were recovered with comparable level for the majority of physiological indexes (Fig 1), except slight difference for POD and SOD activity in leaves (Fig 1D and 1E). Taken together, these results indicated that drought treatment significantly affected the physiology characteristics of M. falcata.

Illumina sequencing and de novo assembly
To further investigate the effect of drought on global transcriptome of M. falcata, total RNAs were extracted from water control plants, treatment plants and recovery plants (with biological triplicates) and subjected to RNA sequencing. In total, we obtained 559,728,990 raw sequencing reads, and 545,584,982 clean reads were obtained after adaptors and low-quality sequences were removed; the value of Q20 is about 97%, Q30 is about 93%, GC content is more than 40.96% (S2 Table). These high-quality reads were then used to assemble transcriptomic data using Trinity. Assembly of clean reads resulted in 475,378 unigenes (S3 Table), ranging from 201-67,238 bp in size, with an N50 length of 1,026 bp and an N90 length of 263 bp (S4 Table). After comparing different transcripts that represent a single unigene, the longest transcript for each unigene was extracted. After eliminating redundant transcripts, a total of 255,014 unigenes were obtained (S3 Table). Among them, transcripts with length longer than 500 bp accounted for about 37.27%, and unigenes longer than 500 bp accounted for about 33.43% (S3 Table). FPKM is an index that were commonly used to evaluate gene expression levels. We  Table), and found that FPKM values of less than 3.57 accounted for a large proportion, but the distribution was unstable in each sample. On the contrary, FPKM values of greater than 3.57 showed a stable distribution in all 9 samples. It is possible that the FPKM value was higher, indicating that the gene expression level was relatively higher, and the expression level is relatively stable in the samples (S5 Table).

Unigenes annotation
Venn diagrams of the number of the non-redundant unigenes in GO, NR, KO, KOG and Swiss-Prot databases were constructed to further annotate transcriptome data (Fig 2, S6 Table). It showed that 147,232 (57.73% of the total) and 113,727 (44.59%) unigenes had significant matches in the Nr database and the Swiss-Prot database, and 204,172 (80.06%) had significant matches in the Nt database. We also found that 113,727 (44.59%) non-redundant unigenes showed similarity to known genes in the Swiss-Prot database. In total, 220,137 unigenes (86.32%) were successfully annotated in at least one of the NR, Swiss-Prot, KO, GO, and COG databases ( Table).

Functional classification of transcriptome data
All transcriptome data were also subjected to three different gene functional classification database, Gene Ontology (GO) database, KOG database and KEGG database. By using Gene Ontology (GO) annotations, total 109,410 non-redundant unigenes were classified into three major functional ontological groups (Fig 3A and S6 Table): biological process, cellular component and molecular function ( Fig 3A). Within the category of biological process, the top three subcategories with the most gene numbers were 'cellular process' (cellular response to stimulus, cell wall organization or biogenesis), 'metabolic process' (hormone, organic substance metabolic process) and 'single-organism process' (single-organism metabolic process, developmental process, cellular process) (Fig 3A). Within the cellular component category, genes involved in 'cell', 'cell part', and 'organelle' were highly represented ( Fig 3A and S8 Table). Additionally, within the molecular function category, genes related with 'binding' and 'catalytic activity' were most highly represented (Fig 3A and S8 Table).
Meanwhile, all unigenes were queried against the KOG database for functional prediction and classification. In total 39,492 non-redundant unigenes were subdivided into 26 classifications (Fig 3B and S6 and S9 Tables). Of these, the 'general function prediction' cluster (including genetic information processing, metabolism, signaling and cellular processes) was the predominant group, followed by 'post-translational modification', 'protein turnover', 'chaperones', 'signal transduction mechanisms', and 'translation ribosomal structure' (Fig 3B and S9  Table).
In addition, analysis on metabolic pathway of the unigenes was conducted using the KEGG annotation system. According to the KEGG database, 50,761 unigenes were classed into 5 KEGG first-level classes (Fig 3C and S6 and S10 Tables), including Cellular Processes (A), Environmental Information Processing (B), Genetic Information Processing (C), Metabolism (D) and Organismal Systems (E). Among them, 19 KEGG second-level classes contained 129 distinct pathways, and the pathways with the largest number of unique transcripts were 'carbohydrate metabolism' (4,288), followed by 'translation' (3,895), and folding, sorting and degradation (3,607) ( Fig 3C and S10 Table).

Analysis on gene expression
In order to discover new candidate genes that may be involved in drought stress in M. falcata, we screened genes that were differentially expressed among control, treatment, and recovery plants. Of 255,014 unigenes, 1,426 differentially expressed genes (DEGs) were identified among control, treatment and recovery group ( Fig 4A, and S4 Table). Among these DEGs, those with higher expression levels in PEG-treated plants than in control plants were considered to be up-regulated (>2 fold), whereas those with lower expression levels in PEG-treatment were considered to be down-regulated (<2 fold). DEGs with higher expression levels in recovery plants than in PEG-treatment plants were also designated as up-regulated DEGs.
In details, 456 and 250 DEGs were up-regulated and down-regulated respectively in treatment vs control, 20 and 35 genes DEGs were up-regulated and down-regulated respectively in recovery vs treatment, and 112 and 615 genes DEGs were up-regulated and down-regulated respectively in recovery vs control (Fig 4B). By clustering genes with the same or similar expression patterns, the function of the unknown gene or the unknown function of the known gene can be identified, therefore, the expression clustering pattern of all DEGs were determined under different experimental conditions (Fig 4C). Many DEGs showed similar expression pattern in both treatment vs control and recovery vs treatment combinations, including genes encoding Appr1 processing enzyme, 30S ribosomal protein S1, trithorax group protein and ADP-ribosylation factor GTPase-activating protein, indicating that drought stress directly induced the expression of these enzyme genes. Meanwhile, we also found that 10 POD genes and one SOD gene were differentially expressed among control, PEG-treatment and recovery group (S3 Fig.), in particular, four POD genes (Cluster91404, 48742, 17838, and 53806) were  Table). There was no significant GO functional enrichment terms in the treatment vs control and recovery vs treatment comparisons. However, in the recovery vs control comparison, we identified enriched DEGs in the cell membrane, suggesting that membrane-related genes preferentially responded to drought tolerance.
Several enriched GO terms related to molecular function were identified, including six binding-related terms (GO:0043169, GO:0046872, GO:0043167, GO:0020037, GO:0046906, and GO:0019825), three enzyme activity-related terms (GO:0016702, GO:0016165, and GO:0003953) in the treatment vs control comparison (S11 Table). We also identified fifteen enzyme activity-related terms, five binding-related terms, and ribosomal structural constituent terms in DEGs from the recovery vs control comparison. With respect to biological processes, enriched GO terms related to development, nodulation, and metabolism were identified in the treatment vs control comparison. Enriched GO terms related to biogenesis, metabolic processes, translation (GO:0006412), cell growth (GO:0016049), carbon utilization (GO:0015976), and gene expression (GO:0010467) were also identified in the recovery vs control comparison (S11 Table). No enriched GO terms were identified in the recovery vs treatment comparison. Taken together, these results suggested that the DEGs in response to drought treatment were more complex during the treatment stage, which included oxidoreductase activity, oxygen-, ion-, and cation binding functions, and oxylipin, alanine, and pyruvate family amino acid metabolism-related genes. Finally, at the recovery stage, DEGs were found to be enriched in oxidase and enzyme inhibitor activity-related GO terms, suggesting that plants showed physiological adaptations to drought.

Transcript factor genes that were induced under PEG treatment
Drought stress also induced the expression of many transcription factor genes in M. falcata (S13 Table). Our results showed that the expression levels of many transcription factor genes, including C2H2 (594) (202) families genes were induced at various level. Among them, fifty-four of them were significant expressed in the treatment vs control, three in the recovery vs treatment, and thirty in the recovery vs control (S14 Table). In addition, expression levels of 54 TFs in treatment vs control were represented by FPKM value (Fig 6 and S14 Table), including AP2, MYB, WRKY, NAC. The AP2-related TF families was the largest families responded to drought stresses identified in M. falcata, which accounted for about one-third of the total number, indicating that AP2 family is directly related to drought stress response.

Drought affects phenylpropanoid and glycolysis/gluconeogenesis pathway
When comparing the gene expression level between the treatment and control groups, we found many genes in the glycolysis/gluconeogenesis pathway were affected. Among them, genes encoding alcohol dehydrogenase class-P (ADH1), pyruvate kinase (PK, phosphoglucomutase), glyceraldehyde-3-phosphate dehydrogenase (gapN), pyruvate decarboxylase (PDC) were up-regulated, and hexokinase (HK) and S-(hydroxymethyl) glutathione dehydrogenase/ alcohol dehydrogenase (frmA, ADH5, adhC) were down-regulated (Fig 7A and S15 Table), which may result in abundant production of ATP and NADPH.

Validation of differentially expressed genes by qPCR
To further validate the transcriptome data, we selected the top 20 genes (up-regulated and down-regulated) with significant differences in each treatment for cluster analysis (S2 Fig). In the treatment vs control comparison, genes encoding NCR peptides, EHBP1 proteins, E3 ubiquitin-proteins, and early nodulin proteins were significantly up-regulated, while genes encoding Na + /H + antiporters, auxin-induced proline rich proteins, ferredoxin-NADP reductase, calmodulin binding protein and ethylene-responsive transcription factors were significantly down-regulated (Fig 5A). In the recovery vs treatment comparison, genes encoding cell wall proteins and NADPH oxidases were significantly up-regulated, while genes encoding ethylene-responsive transcription factors and resistance-like proteins were down-regulated. These results indicated that many genes encoding proteins and transcription factors may function together to coordinate drought stress response.
Furthermore, another 10 genes related to drought stress were also selected for verification by qPCR (Fig 8 and S16 Table). Among them, 5 genes were co-expressed as shown in the venn  diagram ( Fig 4A). The other 5 genes were randomly selected from those that were significantly regulated in the treatment and control combinations, which included genes encoding Myb/ SANT-like DNA-binding domain protein, tocopherol cyclase, zinc-binding alcohol dehydrogenase family protein, auxin induced proline rich protein and ethylene-responsive transcription factor ERF110. Among these ten genes that were analyzed by qPCR, seven of them were verified to have the same expression pattern as from the transcriptome data. In particular, the expression levels of genes encoding trithorax group protein, Appr1 processing enzyme family protein, and 30S ribosomal protein were decreased by PEG treatment but increased after recovery (Fig 8). In addition, genes encoding MYB and zinc-binding proteins were up-regulated by PEG treatment, whereas auxin-related and AP2 genes were down-regulated by PEG treatment (Fig 8). These results verified the accuracy and reliability of the transcriptome sequencing data, suggesting that these genes may function in drought response and they may be explored for alfalfa breeding with drought-resistance.

Discussion
In the present study, we have provided detailed information on the physiological and transcriptome data of M. falcata under drought stress. It can be inferred that M. falcata exhibit stress tolerance mechanism under simulated drought stress. The increase in SOD, POD, and CAT activities (Fig 1), and conductivity and soluble sugar content after stress treatment, indicated that M. falcata possessed a better ROS scavenging ability to resist drought. It was shown that the increase in activity of these enzymes could eliminate ROS and peroxides induced by stress and in turn inhibit the peroxidation of the plasma membrane, and protect cells from damage [38][39][40]. Our results suggest that M. falcata regulates the activity of defense enzymes to reduce the accumulation of harmful substances.
In this study, we also performed transcriptome sequencing and identified 255,014 unigenes in M. falcata (S2 Table). The N50 length of the unigenes was 1,026 bp, and the average length was 654 bp. These results were comparable to those obtained in recently published transcriptomic analyses of other plant species, such as Reaumuria soongorica (N50 = 1,109 bp, average length = 677 bp) [41] and Litchi (N50 = 811 bp, average length = 601 bp) [42]. In addition, more than 80% of the identified unigenes (220,137, 86.32% of the total) were successfully

PLOS ONE
annotated using BLAST searches with the public Nr, Nt, Swiss-Prot, GO, COG, and KEGG databases. This is a considerable proportion, given the absence of genomic information on M. falcata (S7 Table). In the NR database, 75.2% of the gene function was annotated. The most common KOG annotation was "general function prediction only", and this result was consistent with previous studies on alfalfa [43], indicating the great potential for gene function exploration in this plant species. Among these seven databases, 86.32% of the sequences were successfully annotated in at least one database. This means that 13.68% of the obtained gene sequences were not annotated, suggesting that these sequences may be associated with unique physiological processes or stress resistance mechanism in M. falcata (S7 Table). However, incomplete annotation may also result from technical limitations such as low sequencing depth or short read length as reported in other studies [44].
In this study, high-throughput sequencing data was used to characterize the transcriptome of M. falcata, the genomic data of which is currently unavailable. Our study provided the first fully characterized physiological, transcriptome and a comprehensive foundation for future molecular studies on M. falcata. In many studies, digital gene expression sequences could be mapped to the assembled transcriptome for further gene expression analysis. However, transcriptome sequencing results should be verified by qPCR, as was performed in a recent study for Oxytropis ochrocephala Bunge [45]. Many candidate genes involved in drought stress were identified from the RNA-seq data, and we also verified differential expression of 7 specific genes related to drought stress using qPCR. It has long been known that extensive changes in gene expression, including both up-regulated (456 genes) and down-regulated (250 genes) key genes (Fig 4B), occur when plants are exposed to drought stress as in other studies [46,47]. These results suggested that plants varied in their abilities to adapt to drought stress. In extreme cases, gene expression in desert plants may differ extensively from those in inland plants, and differences in gene expression patterns may confer enhanced ability for desert plants in response to drought stress.
In addition, drought-related genes are highly expressed, possibly through increased binding of specific transcription factors (ERF/AP2 and MYB) [48][49][50][51] (Fig 6) to improve drought tolerance, and this finding agrees with previously published works showing that the over-expression of AP2 transcription factor NtERF172 conferred drought resistance [51]. Moreover, alcohol dehydrogenase involved in catalytic decomposition of ethanol, which is eventually oxidized into CO 2 and H 2 O, providing an intermediate product for plant metabolism, eliminate the toxic effects of aldehydes [52], and reduced the toxicity of aldehydes under drought stress [53]. In our study, the gene encoding alcohol dehydrogenase was up-regulated, therefore, it could be the key gene in response to drought stress in M. falcata.
Moreover, GO and KEGG enrichment analyses on DEGs revealed that a number of differential expressed genes were most strongly enriched in biological processes, including citrate cycle (TCA cycle), fatty acid metabolism, carbon fixation in photo synthetic organisms and ribosome, enzymes distributed with a large proportion (S9 Table and Fig 5), indicating a variety of metabolites were synthesized via active metabolic processes in leaves of M. falcata. It has been demonstrated that powerful capability of RNA-Seq technique in identification of novel genes from non-model organisms, and these annotations also provided a valuable resource to study specific processes, and pathways involved in various plant growth and development as in other plant species [54].
Significantly enriched pathways identified by KEGG enrichment analysis include the phenylpropanoid biosynthesis, glycolysis/gluconeogenesis and flavonoid biosynthesis (S9 Table). These important pathways were also found to be significantly regulated under drought resistance in rice [55,56]. Our results also showed that drought stress affected plant respiration, enzyme activity, and metabolic-related processes. It was also revealed that phenylpropanoid biosynthesis was enriched during drought treatment in our study, which was consistent with another research [57,58]. Actually, phenylpropanoids with the greatest potential to reduce ROS were almost exclusively synthesized in response to a plethora of drought stresses [59]. In other studies, plants have resisted external pressures caused by drought stress by increasing antioxidant enzyme (SOD, POD, CAT) activity [60], POD is an essential antioxidant enzyme for plants and it acts synergistically with SOD, CAT, glutathione peroxidase, etc. in a variety of crops to remove harmful free radicals under adverse conditions. Interestingly, six enzymes annotated as peroxidase (POD) of the phenylpropanoid biosynthesis pathway (P<0.05) were up-regulated in treatment when compared to control (Fig 7B). In particular, the increased tolerance to drought stress is related to the enhanced phenylpropanoid pathway, which thereby could increase activity of antioxidant enzymes and reduced membrane damage [8]. Our results demonstrated that M. falcata also possesses a ROS scavenging system that responds to drought stress. The changes of POD reported here were consistent with an increased concentration of antioxidant enzymes, enabling plants to respond to stress and further to adapt to changes in external environment, and this is consistent with increased level of POD ( Fig 1D).
The glycolysis/gluconeogenesis pathway mainly affected ATP supply to respond to abiotic stress [61]. In our study, pyruvate kinase was significantly up-regulated in the treatment and control groups, which is one of the main rate-limiting enzymes in glycolysis (Fig 7A). During glycolysis, pyruvate kinase can convert ADP into ATP, it is presumed that increase in expression level of PK gene may have produced more energy via the glycolysis pathway to tolerate drought [13]. It was reported that accumulation of pyruvate in glycolysis pathway has great impacts, and pyruvate decarboxylase (PDC) controls the anaerobic fermentation of pyruvate. In our study, expression level of PDC gene was significantly up-regulated under drought, therefore, increase of PDC in M. falcata under drought stress not only ensure continuity of the glycolysis process, but also consumes NADH produced during glycolysis process, to avoid cells be damaged by acidification. In summary, efficient up-regulation of genes of entire anaerobic metabolic machinery is essential to provide tolerance against drought [13].
In recent years, increasing research in molecular mechanism were focused on plant responses to drought stress, including studies on gene expression and signaling pathways. It becomes a critical breeding objective to search for drought resistance related genes from M. falcata and apply them in the improvement of alfalfa productivity in many drought-prone environment; it is of great theoretical and practical significance to dissect drought tolerance mechanism in M. falcata. In alfalfa, various regulatory mechanisms were developed to deal with drought, and to avoid or minimize cellular damage caused by drought stress, including higher proline, soluble sugar, malondialdehyde, SOD, POD, CAT content accumulation [38,62]. Furthermore, sequencing strategy was also used to identify drought-responsive micro-RNAs in alfalfa root and leaf tissue [63]. By analyzing the response of M. falcata to drought stress during different stress time, we have detected two significantly enriched pathways, together with 5 genes that are associated with drought resistance. It was evident that phenylpropanoids and glycolysis may play a role in countering drought-induced oxidative damage. This function may be of great value when key components of the highly coordinated network of antioxidant defense are impaired. In addition, the candidate genes can be used as a basis in molecular breeding to improve drought resistance and yields of alfalfa.

Conclusion
In conclusion, we reported the physiological and transcriptome changes in M. falcata under drought treatment and recovery. First of all, our results revealed the relationship between physiological and transcriptome changes under drought stress in M. falcata. Besides, gene enrichment analyses on differentially expressed genes revealed that phenylpropanoid biosynthesis (ko00940) and glycolysis/gluconeogenesis (ko00010) were the two key pathways involved in drought response in M. falcata. Finally, the expression levels of seven DEGs were verified to be involved in drought response in M. falcata by qPCR analysis. Taken together, these results will provide valuable information related to drought stress response in M. falcata and lay a foundation for molecular studies and genetic breeding of legume crops in future research.